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Abstract: We formulate a model of Nf = 4 flavors of relativistic fermion in 2-\-ld in the 
presence of a chemical potential fi coupled to two flavor doublets with opposite sign, akin 
to isopsin chemical potential in QCD. This is argued to be an effective theory for low 
energy electronic excitations in bilayer graphene, in which an applied voltage between 
the layers ensures equal populations of particles on one layer and holes on the other. The 
model is then reformulated on a spacetime lattice using staggered fermions, and in the 
absence of a sign problem, simulated using an orthodox hybrid Monte Carlo algorithm. 
With the coupling strength chosen to be close to a quantum critical point believed to 
exist for Nf < Nfc ~ 4.8, it is found that there is a region below saturation where both 
the carrier density and a particle-hole "excitonic" condensate scale anomalously with 
increasing /i, much more rapidly that the corresponding quantities in free field theory, 
while the conventional chiral condensate is strongly suppressed. The corresponding 
ground state is speculated to be a strongly-correlated degenerate fermion system, with a 
remnant Fermi surface distorted by a superfluid excitonic condensate. The model thus 
shows qualitatively different behaviour to any model with /i 7^ previously studied by 
lattice simulation. 
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1 Introduction 



The study of quantum systems containing a non-zero density of conserved charge be- 
yond perturbation theory remains a challenging but fascinating problem. Estimating 
Euclidean Green functions via Monte Carlo importance sampling of the discretised path 
integral is rendered inoperable once /x/T ^ 1, where fi is the chemical potential asso- 
ciated with the conserved charge such that the charge density n = d\nZ/dfi, because 
generically the integrand e"'^-^*-'^'', with Se the action in Euclidean metric, is not positive 
definite once fi 0. There remain a small number of interesting problems with real Se 
which can be tackled by orthodox methods: fermionic models such as Gross-Neveu and 
NJL, which with fi are either Fermi liquids [Ij or weakly-interacting BCS super- 
fluids [21 [3]; certain gauge theories with groups such as SU(2) |3] or G2 [5] containing 
real matter representations; and QCD with non zero isospin density [6]. There is also a 
recent study of the interacting relativistic Bose gas which deals with a complex Se by 
sampling a complexified configuration space using Langevin dynamics [7j. 

In the present paper we claim to add to this list a model of strongly-interacting 
fermions in 2+ld inspired by recent developments in graphene. Recall that due to 
the special properties of the honeycomb lattice (with just two isolated zeros or "Dirac 
points" of the dispersion E{k) in the first Brillouin zone), low energy excitations of either 
electron or hole nature are naturally described by a 2+ld Dirac equation. For monolayer 
graphene, the counting of degrees of freedom (2 Dirac points x 2 C atoms per unit cell x 
2 electron spin components) results in Nf = 2 flavors of 4-component relativistic spinor 
fields. Here we consider a model of bilayer graphene comprising Nf = 4 relativistic 
flavors, in which a biassing voltage is applied across the layers so that a non-zero density 
of electrons is induced on the negative sheet and an equal density of holes induced on the 
positive. It will be shown in Sec. [2] that this is equivalent to introducing an equal and 
opposite chemical potential fi on each sheet, so that the resulting relativistic effective 
theory has an "isospin" chemical potential fij with two "m" and two "d" flavors. Just 
as for QCD, it is straightforward to show the resulting Se is real and positive so that 
Monte Carlo simulation is applicable. 

The role for numerical simulation becomes apparent once inter-electron interactions 
are considered. The Coulomb interaction between charges in graphene is effectively en- 
hanced by a factor vf/c ^ 300, where vp = dE / dk\k=kp is the Fermi velocity of the 
charge carriers, which is constant for a linear dispersion. Since in undoped graphene 
with carrier density = only quantum fluctuation (namely vacuum polarisation) ef- 
fects contribute to screening, the resulting effective theory of relativistic electrons has 
an effective fine structure constant ~ 0(1), albeit with a non-covariant "instanta- 
neous" interaction of the form AqiP'^qiP due to vp ^ c. This implies that the theory 
must be treated as strongly-interacting, and the possibility of non-perturbative effects 
such as particle-hole condensation {ipip) 7^ 0, resembling chiral symmetry breaking and 
resulting in djTiamical mass gap generation leading to a Mott insulator phase, should 
be considered. Son [8], in the context of a model with Nf relativistic species and vari- 
able coupling strength a has suggested that gap generation occurs for Nf sufficiently 
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small and a suffciently large. Moreover, the critical coupling a^Nf) defines a quantum 
critical point (QCP) where the scaling of operators and correlation functions in princi- 
ple receive anomalous corrections, which could for instance modify the linearity of the 
electron dispersion. 

The scenario was tested in a Monte Carlo simulation with 2+ld staggered fermions 
(at a QCP the underlying lattice should become irrelevant) [9j; in the strong coupling 
limit a — 7- oo the QCP is found for Nj^ = 4.8(2). Simulations incorporating a more 
realistic long-ranged Coulomb interaction estimated ac = 1.11(6) [10] for the case Nf = 2 
relevant for monolayer graphene, suggesting there is a real possibility of gap generation 
for suspended samples (the effective value of a is sensitive to the dielectric properties 
of the substrate, and is expected to be maximal when the substrate is absent). To 
date there is no definitive experimental signal for chiral symmetry breaking or mass gap 
generation; however, non-linearities in the electron dispersion as a result of interactions 
have been reported in [II]. Even if gap generation does not take place, it is possible 
that physical graphene lies close to the QCP in the chirally symmetric phase, implying 
modifications to electron transport. 

For the case of a voltage-biased bilayer a new condensation channel, between particles 
in one layer and holes in the other, opens up; in what follows we will refer to this as 
"excitonic" condensation. Because of the increasing density of states, as the Fermi energy 
is raised this should become the preferred channel even though inter-layer interactions are 
weaker than intra-layer ones. Once again, gap formation results implying an insulating 
phase; since this can be controlled by the external voltage the possibility of graphene- 
based electronic components arises [12] . A self-consistent treatment however, taking into 
account the enhanced screening due to Nf = 4, finds the resulting gap A//i ~ O(10~^), 
suggesting that excitonic condensation will be difficult to achieve experimentally [13]. 

Our purpose in this paper is to re-examine excitonic condensation in the non-perturbative 
setting required by the vicinity of a QCP. We will reformulate the lattice model used for 
the investigation of the QCP in undoped graphene in [9l HI] to include Nf = 4 contin- 
uum Dirac fermions and a nonzero isospin chemical potential /x to represent the voltage 
bias. In order to formulate the model with a local interaction on a 2+ld lattice, the 
interaction between fermions is forced to have the form of a contact between local charge 
densities, schematically {ip'yoipy, and hence no long-range Coulombic tail. As argued in 
[9], we expect that large vacuum polarisation effects make this irrelevant near the QCP 
(or equivalently in the large Nf limit). The model is presented in detail in Sec. [2l An 
important and necessary simplification is that inter- and intra-layer interactions between 
fermions have a common coupling and hence the same stength. Since n ^ boosts the 
density of available particle-hole states, we again expect this to become less relevant as 
fi is increased. 

In Sec. [3]we present results from numerical simulation of the model on lattice volumes 
32^ and 48^. Since it supposedly describes a continuum effective theory in the vicinity 
of a QCP it is diffcult at this exploratory stage to assign physical units to the simulation 
results, or even in principle to know the physical anisotropy ratio at/ as- We choose a 
coupling strength close to the putative QCP for Nf = 4; since the strong-coupling limit 
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is hard to isolate for our formulation [15] this proves to be a non-trivial exercise, and 
indeed we will see it appears our simulations lie just in the chirally symmetric phase at 
fi = 0. Next, we introduce 7^ and monitor the chiral condensate, carrier density and 
exciton condensate as it is increased. It will be shown that excitonic condensation does 
indeed take place and that both carrier density and exciton condensate considerably 
exceed the values expected for free fermions (augmented by a small symmetry-breaking 
source term), whereas the chiral condensate is suppressed. We discuss our findings in 
Sec. m Although the model is motivated by condensed matter physics, we will argue 
the results are of wider interest, and yield perhaps the first insight into Fermi surface 
physics in the presence of genuinely strong interactions. 



2 Formulation and Interpretation of the Model 

Here we outline the formulation of an effective field theory for the graphene bilayer. 
Physically, the idea is that there are Nf = 2 flavors of relativistic fermion on each 
monolayer, described by an action in Euclidean metric fl6[ [8]: 

'S'mono = ^ I dXoSx{ll)alodollJa + Vplpa^-Vlpa + i^oV^aToV'a) + ^ / dXod^ x{diAof , 

(1) 

where e is the effective electron charge (whose value depends on the dielectric properties 
of the substrate), and the 4x4 Dirac matrices satisfy {7^, 7,^} = 25^^, = 0, 1, 2. Note 
that for this reducible represention of the Dirac algrebra there exist two independent 
matrices 73 and 75 = 7o7i7273 which anticommute with the 7-matrices present in ([1]). 
Aq is a fluctuating 3 + Irf electrostatic potential field sourced by the charge density i''jo4', 
and is a remnant of the full electromagnetic field in the instantaneous approximation 
justified for vp <^ c. 

Now, for a perfect bilayer formed from two monolayers stacked in AB configuration 
with interlayer coupling strength t' ~ 0(0. l)t where t is the hopping parameter in 
the monolayer tight-binding Hamiltonian, it is known that the dispersion relation for 
massless fermions in the low-energy limit in the vicinity of the Dirac point is quadratic, 
only becoming approximately relativistic (i.e. linear) for ka ^ t'/t \i7\. For fi ^ t 
the dispersion then takes the expected form = (yU ± vpk)'^ [H]. However, recent 
theoretical studies of strained bilayers suggest that under mechanical deformation the 
parabolic bands split to form separate Dirac cones, so that in this case a description in 
terms of Nf = 4 relativistic species is not a bad approximation [19] . Our formulation 
makes the additional, perhaps unwarranted, approximation that interactions between 
charge carriers on different layers are of identical strength and character to interactions 
within a layer - the necessity for this will become clear below. 

The second ingredient of the model is that the layers are given equal and opposite 
constant bias voltages ±/i, inducing on one layer a negatively charged concentration of 
particles and on the other a positively-charged concentration of holes. As the notation 
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implies, the bias voltage is equivalent to a chemical potential, and in fact the theory is 
formally very similar to the case of QCD with isospin chemical potential fij = fii = — /i2, 
where the subscripts which here label the layers usually stand for the light quark flavors 
u and d. Euclidean formulations of systems with fi are generally afflicted with a 
"Sign Problem" , ie. the Lagrangian density C is no longer positive definite, or even real, 
since the inequivalence under time reversal translates into inequivalence under complex 
conjugation in Euclidean metric. This makes Monte Carlo importance sampling as 
a means to handle strongly fiuctuating observables inoperable. However, the case of 
isospin chemical potential is known not to have a Sign Problem and is hence simulable 
using orthodox methods, as we shall now demonstrate. 

If we denote the fermion degrees of freedom on one layer by ip and on the other by 
0, define units so that Vp = 1, and write ^^=q i 2 9^"^^ + [iAq + /i)7o = D[A; ji], then 
the fermion part of the Lagrangian can be written: 

(2) 

Here we have introduced two new real parameters: m is an artificial bare mass which 
induces a gap in the fermion dispersion relations, and whose sign has no physical conse- 
quence for a single fiavor in the absence of interactions; j a source strength coupling ip to 
0, thus linking the layers and eventually enabling calculation of the exciton condensate. 
In principle both m — )■ and j — )■ limits need to be taken in order to make contact 
with physical bilayer graphene. Integration over the Grassmann bispinors ^ then 
results in the functional measure detA^fA]. 

An important identity which the model inherits from the gauge theory is 

D^[A-^i] = -D[A--^i]. (3) 

It is then straightforward to check (assuming the dimension of D is even) that 

detM=dei[{D + m){D + m)^ +f]>Q, (4) 

and 

\ [D + Tn){D + my + j'^ J ^ ' 

implying both that 

^QiM^M = det^M, (6) 

and also that the desired functional measure det results from integrating over bosonic 
fields $ starting from a non-local "pseudofermion" Lagrangian 

Cpf = <^^[{D + my{D + m)+f]-'<^. (7) 

This has the practical advantage that $ has half as many degrees of freedom as and 
makes eqn. ([7]) the appropriate starting point for the design of a hybrid Monte Carlo 
simulation algorithm. 
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The specific version of D+m in our lattice model employs single-component staggered 
fermion fields ipx, 4>x defined on the sites of a 2 + Irf square lattice, with a non- compact 
formulation of the electrostatic potential formally defined on the link joining sites x 
and X + 0: 

i=l,2 

(8) 

where the signs r^^a; = (— l)^""* ^^^-i ensure Lorentz covariance in the long wavelength 
limit. It can be shown that the relation between the number of staggered fields N 
(counting ilj,(f) yields N = 2) and the number Nf of continuum Dirac 4-spinors is |20] 

Nf = 2N. (9) 

It is worth noting the global symmetries present in the model. For fi = m = j = 
the continuum action (j2]) is invariant under a U(8) rotation \E' H- U"^, \l/ i— ?• where 

= 2^^7375. This symmetry is broken to (U(4))^ by /i 7^ 0, and then to (U(2))^ by 
m 7^ 0. Setting the interlayer coupling j ^ with m = locks the ip and components 
together, so that in this case the residual symmetry is U(4). For the staggered lattice 
fermions of ([H]) the original symmetry is U(2)®U(2)£, where the second rotation can be 
written as U{a] x) = exp{iexC(°'A°'), where A" is one of the four hermitian generators of 
U(2) and = (-l)^o+xi+x2_ getting ^ breaks this to (U(l)(g)U(l) followed by 
m ^ 0,j = to (U(l))2, and m = 0,j ^ to U(l)(g)U(l)^. 

The fermion action is supplemented by a Gaussian weight for the A fields: 

where is a parameter governing the strength of the coupling between the potential and 
the fermions. The resulting dynamics describes A fluctuations having the same form as 
the continuum action ([T]) in the strong-coupling or large- A^j limits, but for which explicit 
screening removes the long-ranged tail away from these limits; further justification 
for this approximation is given in [9l[Tl]. For Nf = 2 this formulation yields an identical 
path integral to the lattice action couched in terms of compact link variables given in 
Eqn. (7) of [11]. For Nf > 2, however, the two approaches are not equivalent since the 
compact formulation leads to extra terms of the form {ipipclxpY in the effective action - 
although these operators may well be irrelevant at the critical point. The exact lattice 
version of the non-compact action for /i = and arbitrary Nf once A is integrated out 
is given in eqn. (2.2) of [2T] . 

We next discuss the implications of relaxing the requirement that inter- and intralayer 
interactions between fermions are identical. The non-trival terms in the action are of 
the form tpUe'^tp, il)U*e~^il), where f/ is a complex number not constrained to have unit 
modulus. Integration over U leads to repulsive particle-particle and hole-hole interac- 
tions, and attractive particle-hole interactions. Suppose we wanted to make the model 
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more realistic by introducing a distinction between intralayer and interlayer interactions. 
One way to do this would be to introduce a second boson field coupling to ip and with 
opposite signs, in effect introducing repulsion between ■j/'-particles and 0-holes so that 
the ip-cf) and (p-ip couplings are weaker than those of ip-ip or (p-cj). The interaction terms 
could then be written ^/^^f/l^e'^^^+o, (t)xUV*e-y^^(^, -^^U*V*e-^'^lj^_Q, -(l)^U*Ve''(j)^_o, 
etc. In the limit V — > 1 integration over ip, ip leads to a factor detD[fi\, while integration 
over 0, (f) gives detD[—fi\. With the help of (|3]) we confirm the resulting functional mea- 
sure detD[/i]D^[/i] is positive definite. In the limit f/ — )■ 1, however, the same process 
leads to detD[^]D*[—^] = det^Z)[/i], which is no longer positive definite. In other words, 
attempting to make the model more realistic reintroduces a Sign Problem, although a 
more detailed study would be needed to determine its severity. 

Now let's discuss observables. The usual chiral condensate (which has been called 
the "exciton condensate" in our earlier work [9l [H] ) is given by 

(M) = ^ = (#)-(#). (11) 
om 

Note the sign of the condensate is not physical, and that the two terms on the RHS of 
( ITTl) give equal contributions. From the discussion above it should be clear that for /i 7^ 
formation of this condensate spontaneously breaks (\J{l)<S)U{l)eY to (U(l))^, resulting 
in two Goldstone modes in the limit m — )■ 0, j — )■ 0. The exciton condensate discussed 
in [T3] and which is the main focus of this paper is given by 

{^^) = ^— = i{^<P-^^). (12) 

In this case the symmetry breaks to U(l)®U(l)e implying the same number of Gold- 
stones. In fact for /i = and m = j, ('I'\E') and (\E'\Ef) are physically indistinguishable, 
both being equivalent to the chiral condensate of the Nf = 2 theory. Figure [1] below 
confirms that with /i = our code generates results consistent with (\E^\E') = y. 

With /X 7^ we next define the charge carrier density 

= = (tpDoij) - {4>Do4>)- (13) 

Once again, both terms on the RHS give equal contributions - the first term represents 
the density of electrons in layer 1, and the second the density of holes in layer 2. 

Figure [1] shows the results of a pilot run on 8^ at = 0.4 and ma = 0.05. For 
ja = 0.05 the two condensates are degenerate at = as argued above. As fi increases, 
our naive expectation is that a Fermi surface of radius fi forms on each layer, one 
containing particles, the other holes, implying ric oc fi"^. As /i grows, ipip and (fxp pairing 
are suppressed because a free particle-hole pair costs energy 2/x to create at either Fermi 
surface, whereas ipcp pairing is promoted, because it costs zero energy to create a particle 
on one Fermi surface and a hole at the other, with the density of states at either increasing 
oc /i. Thus (^f\Ef) decreases as /i rises from 0, while (\E'\E') increases. The rise in (\E'^') 
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Figure 1: Fermion condensates as a function of /i at g — 0.4 on 8'^ with bare mass ma = 0.05, 

ja = 0.05,0.1. 



seems to be relatively more pronounced for smaller j. This trend persists until /xaj ^ 0.3. 
What happens after that should be understood in terms of saturation, an artifact which 
sets in when the fermion density is a significant fraction of its maximal value of one per 
lattice site. With our normalisation of ric this sets in for ^at ~ 0.5, a surprisingly small 
value based on experience with other models. In a saturated world fermion excitations 
of all kinds are kinematically suppressed, and the condensates tend to zero in this limit. 



3 Numerical Results 

Our strategy in this paper is to investigate the effect of varying in our bilayer model 
flHfTOj) starting close to the quantum critical point. The first task is to find the coupling 
g1 where the QCP is located for Nf = 4; we use the approach [151 IS] of searching for a 
maximum of {^'^) as is varied and identifying that with the strong coupling limit 
of the continuum model. We then assume g~'^ ^ S'peak' since if the value Nfc = 4.8(2) 
obtained in [9] is universal there should only be a narrow range of g'"^ corresponding to 
the chirally broken phase. The results for (^'^(m)) in Fig. |2]show that g'^^^y- ~ 0.30, 
much larger than the value ~ 0.05 obtained with the compact formulation [S]. Another 
contrast with previous work is that it is also apparent that g~^^-^^ increases with m, from 
roughly 0.275 at ma = 0.07 to 0.35 for ma = 0.01, although at this stage we cannot 
exclude the possibility that finite volume effects influence the result. For small m a 
linear extrapolation to the chiral limit seems reasonable; we conclude, conservatively, 
that in this limit g~^^^ e (0.275,0.35). 

Figure [3] shows {^"^) data as a function of m for g~'^ fl'peak- Whilst the quadratic 
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Figure 2: (Color online) (5'^') vs g ^ for iV/ = 4 and various m near gpjfak ~ 0.30. The simulations 

were performed on both 32'^ and 48^^ lattices. 
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Figure 3: (Color online) {^'^) vs m for g ^ = 0.35, 0.375, 0.40 fitted to a quadratic polynomial. 



extrapolation to the chiral limit is not conclusive, the marked non-linearity of the fits 
suggests the QCP value g'"^ lies close to this region; however, a much more extensive 
simulation campaign would be needed to pin it down. For our purposes it suffices to work 
close to a strongly-interacting QCP, while leaving the issue of whether chiral symmetry 
spontaneously breaks unresolved. Henceforth, all numerical results are obtained with 
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varied. Unless otherwise stated, the chiral limit m = will be assumed. 
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Figure 4: (Color online) (W^E'} vs ^ on 32^ for m = and ja — 0.01,0.02,0.03. Dashed lines show the 

same quantity evaluated for free fields. 
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Figure 5: (Color online) vs /i on 32^ for ja = 0.02 and ma = 0, 0.01, 0.02, 0.03. The dashed line 

shows the same quantity evaluated for m = for free fields. 

Figure H] shows the exciton condensate (\E'\E') as a function of /i for three different 
j. The figure shows the same broad features as Fig. [H namely a rapid rise to a fairly 
sharp maximum at /xa ~ 0.3, followed by a still more rapid fall; the signal is very small 
indeed by /la = 0.6. As we shall see in Fig. |H1 at this value of /i the system has reached 
saturation with a maximum possible density of particle-hole pairs consistent with the 
Pauli exclusion principle on a fixed lattice; our model can only be interpreted as a 
description of bilayer graphene for values of /i much smaller than this. 
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The dashed hnes in Fig.lHshow (\E'\1/) evaluated using the same measurement code but 
with g"^ set to zero, yielding the value for free fields. Since the (U(l)®U(l)e)^ symmetry 
is manifest for j = the free-field condensate must vanish in this limit, and the curves 
are consistent with this expectation. The large disparity between (\E'\E')int and (\E'\E')free 
notable in the range 0.2 ^ /xa ^ 0.4 signals that (U(1)®U(1)£)^ is surely spontaneously 
broken here. Close inspection of the figure reveals that (\Ef\E')free rises monotonically, 
but not quite smoothly, with /i until reaching a maximum at /la ^ 0.9. The disparity 
with the apparent saturation observed in the interacting model will be further discussed 
below. The barely visible wiggles are probably a finite volume artifact similar to that 
noted in studies of another system with a Fermi surface [3] . Figure |5] plots the same 
scan but this time showing that the effect of varying m is negligible except for the very 
smallest values of /x. Since the operator \E'\E' is constructed to be conjugate to j, not m, 
this is as expected. 
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Figure 6: (Color online) (^'\E') vs j for m = and various ji on 32^ (open) and 48'^ (closed symbols). 
Dotted lines show fits to eqn. (fT4| . Dashed lines show the same quantities evaluated for /i = 0,0.2 on 

48'^ for free fields. 

In order to interpret the condensate data it is necessary to extrapolate j — )■ 0. 
Figure [6] shows for several j on two different volumes, together with extrapolations 
of the form 

{mm) = {-^-^(j = 0)) + Aj + Bf + Cf. (14) 

Taking finite volume effects into account, it seems that at least for fia > 0.10 the 
fitted intercept is non-vanishing, confirming the spontaneous breaking of particle-hole 
symmetry due to excitonic condensation (\E'\E') ^ 0. The extrapolated condensate is 
shown fitted to a power law of the form (\E'\E'(j = 0)) = ai/i"^ in Fig. [3 the fitted 
parameters are 

ai = 7.0(2); aa = 2.39(2). (15) 

The power-law rise is more rapid than would be expected from a BCS-style mechanism 
driven by condensation of particle-hole pairs in the immediate vicinity of a Fermi surface. 
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Figure 7: (Color online) (^^'(j — 0)) vs /i on 48^ fitted to a power law for fi = 0.05 — 0.20. The dashed 

line corresponds to exponent 02 = 2.39(2). 



This is because in a BCS condensation the density of available pairing states scales with 
the area of the Fermi surface, oc yU*^"^ in d space dimensions. Despite this somewhat 
empirical approach, the non-linear increase of with /i is a robust conclusion at 

variance with a conventional weakly-interacting BCS scenario. 




Figure 8: (Color online) Carrier density ric vs /i on 32'^ 



and j = 0.01,0.02,0.03. The dashed 



line shows the same quantity evaluated with j — 0.01 for free fields. 



Next we consider the carrier density Uc defined in ( fT3|) . and shown in Fig. [8l This 
rises monotonically from zero with n until /za ~ 0.5, when saturation sets in; the effect 
of j 7^ is to round off this behaviour by reducing the carrier susceptibility \dnc/dfj,\ 
slightly. Once again, the contrast with the free-field behaviour, which only reaches 
saturation at yua ~ 1.3 and is shown by the dashed line, is marked. 
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How should we interpret the finding that n™* ^ n^^^^l For degenerate fermions 
the carrier density, remembering to count both particle and hole states, is given by 
Tic = kp/2'n'. For free massless fermions the Fermi energy /i is equal to Fermi momentum 
kp] If we wish to retain the notion of a Fermi surface (albeit one distorted by exciton 
condensation) for the interacting system, we are forced to conclude n ~ Ep < kp imply- 
ing strong self-binding, i.e. the degenerate fermions have a large negative contribution 
to their bulk energy. This is a feature of working near a QCP, and was not observed, 
e.g. in studies of relatively weakly-interacting systems at non-zero density such as the 
Gross-Neveu model in 2+ld [22j where interactions are suppressed by '^/Nf, or two color 
QCD |1] where the quark density Uq ^ n^^^^ all the way to saturation. 
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Figure 9: (Color online) Carrier density ric vs j on 48'^ for various /i. Dotted lines show a quadratic 
extrapolation j 0. Dashed lines show the same quantity evaluated for free fields with ^ = 0.1,0.2. 

As before, the region of physical interest is for /i well below saturation: Figure [9] 
plots the variation of with source strength j, togther with a quadratic extrapolation 
to j = 0, showing that the effect of j 7^ for this observable is regular but certainly not 
negligible. Finally Fig. [TU] plots nc(/u; j = 0) together with a power law fit = 
The fitted parameters are 



bi = 18.6(4); 62 = 3.32(1) 



(16) 



As expected, the fitted value of 62 considerably exceeds the naive expectation Uc oc /i^ 
based on a weakly-interacting system. 

In Fig. [TT]we show the chiral condensate order parameter as a function of fi for 

various values of the symmetry breaking parameter m at fixed ja = 0.02. Its magnitude 
at /i = falls approximately linearly with m implying restoration of chiral symmetry as 
m — 0. Even so, it exceeds the free-field value by over a factor of two, refiecting the 
vicinity of the QCP. Note though that a measure of the density of particle-hole 
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Figure 10: (Color online) Carrier density ndj = 0) vs /i on 48^ fitted to a power law for [i = 0.05 — 0.20. 
The dashed line corresponds to exponent 62 = 3.32(1). 
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Figure 11: (Color online) Chiral condensate vs [i on 32^ for j = 0.02 and m = 0.01,0.02,0.03. 

Dashed lines show the same quantity evaluated for free fields. 



pairs in the condensate, is roughly one-third of the peak value of the exciton condensate 
|(\E'\E')| seen in Fig. HI As fi increases K^^^P)] falls monotonically reflecting the increasing 
difficulty of particle-hole pairing within a layer as the biassing voltage rises. We deduce 
that near the QCP the impact of the biassing voltage is to favour inter-layer over intra- 
layer pairing; indeed the inter-layer pairing is suppressed completely, falling below even 
the free-field value, by fia = 0.5 where saturation sets in. 

Finally we report on a preliminary calculation of the spectrum of quasiparticle exci- 
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Figure 12: (Color online) Normal and anomalous fermion energies E{k — 0) ys j on 48'^ for ^, 

0,0.1,0.2. 



tations, obtained from analysis of the following fermion correlators: 

C^ik.t) = ^(V^(0,0)V^(f,t))e- 

X 

CA{k,t) = ^(^(O,O)0(x,t))e- 



-ik.x. 



ik.x 



(17) 



Due to the form of the staggered fermion action ([H]) the set of two-dimensional sites 
X only includes those displaced from the origin by an even number of lattice spacings 
in any direction, and the physically accessible momenta have ki = 2Tini/Ls with rij = 
0, 1, . . . , Ls/4. We distinguish between the normal propagator Cn describing carrier 
motion within a layer, and anomalous propagator Ca describing inter-layer hopping, 
which relates eg. destruction of an electron on layer 1 to creation of an anti-hole on 
layer 2. On a finite system Ca is non- vanishing only for j 7^ 0. 

In this first study we have considered k = Q only. In accordance with a study of 
quasi-particle propagation in a thin-film BCS superfiuid [2] we find that the correlator 
signal resides in Re(C7v) and Im(CA), and that in the chiral limit m — )■ CAr(t) = for 
t even and Ca(^) = for t odd. We thus fit the correlators on every second timeslice for 
the excitation energy E using the forms 



Re{CN{k,t)) 
lm{CA{Kt)) 



(18) 



where in general A ^ B ioi fi 0. The resulting energies are shown for small /x as 
a function of j in Fig. [121 Two features are apparent: firstly normal and anomalous 
channels yield consistent results, as expected [2], although the normal data have smaller 
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errorbars; secondly the extrapolation j — )■ appears smooth, and suggests limj_j,o E{j) > 
for /i 7^ 0. In other words, a voltage bias induces anomalous propagation indicative of 
particle-hole mixing, a manifestation of excitonic condensation 7^ 0. This should 

also result in a non- vanishing energy gap at the Fermi surface, but confirmation requires 
a study of propagation with k ^ [3\. 

4 Discussion 

In this paper we have set out an effective (albeit simplified) field theory for low-energy 
charge transport in voltage-biased bilayer graphene, and shown how it can be simulated 
using orthodox lattice field theory methods, because its action in Euclidean metric is 
real. An interesting feature of the numerical formulation is that it is possible to run in 
the chiral limit m = so long as the ip(f) coupling j 7^ 0. There are formal similarities to 
QCD with non- zero isopsin density |6]; however the resulting dynamics differs sharply. 
While QCD is an asymptotically-free theory implying that eventually a weakly-coupled 
description becomes valid as /x — )■ 00, here the field correlations remain strong at all 
scales, even in the absence of confinement, due to the vicinity of the QCP. For this 
reason the model is of intrinsic theoretical interest independent of any possible physical 
applications for graphene. 

Precise location of the QCP by numerical means has proved challenging; nevertheless 
the curvature of the chiral condensate data {^"^{m)) of Fig. [3] are suggestive of a critical 
scaling {^^) oc m« expected at or near a QCP. Equation of state fits predict S in the 
range 2.7 (the value for monolayer graphene with Nf = 2) [13] to 5.5 (the value in the 
strong coupling limit Nf = Nfc ~ 4.8(2)) [9]. Considerably more work would be needed 
to confirm this quantitatively. 

Our main results in this first study are therefore qualitative. Runs with j 7^ yield 
measurements of the exciton condensate which show a rapid rise as /i is increased 

from zero. The data extrapolated to j — )■ suggest that condensate remains non- 
vanishing in this limit consistent with spontaneous symmetry breaking and superfluidity; 
indeed the data of Fig. [7] permit a power-law fit (\E'^) oc /i^'^. This is notable because 
a weak-coupling BCS description of superfluidity predicts the condensate should scale 
with the area of the Fermi surface, namely oc /i. Similarly, the carrier density 

ric oc /i^'^ (Fig. in]), to be contrasted with the weak-coupling behaviour oc /i^. With the 
resolution we are working with there is no sign of an onset value of the chemical potential 
/io > 0, such that ric = 0, (\&\&) = for < /Iq. This is another important contrast with 
the systems studied in [H [21 [S] HI E] |6] . The likely reason is that at the couplings studied 
there is no mechanism for spontaneous mass generation, so that the lightest degree of 
freedom carrying a conserved charge remains massless. The final interesting observation, 
shown in Fig. [TTl is that the chiral conndensate is strongly suppressed as /i rises, 

presumably because of the rapidly-increasing energy cost of a particle-hole pair within 
a layer, and is consistent with zero post-saturation. 

Another observation to note is that below saturation both 11^ ^ n^^'^'^ and (\E'\E') ^ 
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(^xl/xl/^frcc Q^(^g again, this is indicative of strong correlations among the fields, such that 
Ep < kp, as is the precocious value of fia at which saturation sets in. It suggests that 
the self-consistent diagrammatic approach of [I3j (which employs large- A^/ methods so 
does not probe the QCP) may yield an unduly small estimate of the condensate. It must 
be stressed, however, that in the absence of a physical scale setting any phenomenogical 
applications of the model to real graphene are premature. 
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Figure 13: (Color online) The ratio vs fia on 32^ (filled) and 48"^ (open) for jc 

0.01,0.02,0.03, together with the j ^ extrapolation on 48^. 



In conclusion we claim to have initiated a lattice Monte Carlo study of strongly 
interacting degenerate fermions, which displays significant qualitative differences to other 
degenerate systems studied previously. A final question worth discussing is to what 
extent the concept of a Fermi surface, either sharp or distorted by particle-hole excitonic 
condensation, remains intact in a strongly- interacting environment. Departures from the 
canonical weak-coupling are manifested as anomalous scaling with Fermi energy fi (see 
eqns. f|T5][T6|) ): however recall that in an interacting Fermi liquid the relation between 
particle density and Fermi momentum kp, namely Uc oc kp, should remain inviolate. In 
the BCS picture, the density of condensed particle-hole pairs (\I'\I') arising from plane 
wave states within a shell of thickness A around the Fermi surface implies 

(^^) oc Akp^ oc AnP^. (19) 

To test whether the scaling ( IT9l) is retained even at strong coupling. Fig. [13] plots the 

ratio (\E'\E') / ric vs. fi for various j on two volumes, together with the j — )■ extrapolation 
on 48^. It looks plausible for fia ^ 0.2, on the assumption that the gap A has a strong 
/x-dependence, which should be the case for near-conformal dynamics. It may well prove 
possible, therefore, to define a Fermi surface in the vicinity of a quantum critical point. 
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